require(readxl)
require(bizdays)
library(foreign)
require(splitstackshape)
#require(haven)

is.date = function(x) inherits(x, 'Date')

create_threat_data = function(this_year, base_file){

	admin_buffer = 0
	me_buffer = 0
	
	d = base_file[which(base_file$year == as.character(this_year)),]
	d2 = school_base_file[which(school_base_file$year == this_year),]
	
	academic_calendar = readxl::read_xlsx(path=paste(input_gdrive,"calendar_", this_year,".xlsx",sep=''),sheet = "Sheet1")
	
	holidays = as.Date(academic_calendar$holidays[!is.na(academic_calendar$holidays)])
	important_dates = academic_calendar[!is.na(academic_calendar$important),c("name","important")]
	important_dates$important = as.Date(important_dates$important)

	cal = create.calendar(name="evaldays", holidays=holidays, weekdays = c("saturday","sunday"), 
		start.date = important_dates$important[which(important_dates$name == "First Day")],
		end.date = important_dates$important[which(important_dates$name == "Last Day")])

	cutoff_m1 = important_dates$important[which(important_dates$name == "Master Educator Cutoff 1")]+me_buffer
	cutoff_p1 = important_dates$important[which(important_dates$name == "Administrative Cutoff 1")]+admin_buffer
	cutoff_p2 = important_dates$important[which(important_dates$name == "Administrative Cutoff 2")]+admin_buffer
	eval_end_date = important_dates$important[which(important_dates$name == "Evaluation End Date")]
	dccas_end_date = important_dates$important[which(important_dates$name == "DCCAS Testing End")]
	dccas_start_date = important_dates$important[which(important_dates$name == "DCCAS Testing Start")]
	eval_start = important_dates$important[which(important_dates$name == "Evaluation Start Date")]

	end_winter_break = academic_calendar$holidays[which(academic_calendar$name == "Winter Break End")]

	bizdays(important_dates$important[which(important_dates$name == "Evaluation Start Date")], dccas_start_date, cal)
	bizdays(cutoff_m1, dccas_end_date, cal)

	#################################
	#
	# ME NO THREAT - 2 CYCLES
	#
	#################################

	# ME cycle 1
	late_tests = which(d$obsdate_1m > cutoff_m1)
	nothreat_m1 = rep(NA, nrow(d))
	nothreat_m1[-late_tests] = bizdays(d$obsdate_1m[-late_tests], cutoff_m1, cal)

	max_ntm1 = rep(bizdays(eval_start, cutoff_m1, cal), nrow(d))

	max_ntm1_postbreak = bizdays(end_winter_break, cutoff_m1, cal)
	ntm1_prebreak = nothreat_m1 -  max_ntm1_postbreak
	ntm1_prebreak[which(ntm1_prebreak < 0)] = 0
	ntm1_postbreak = nothreat_m1 - ntm1_prebreak
	
	# ME cycle 2
	nothreat_m2 = rep(NA, nrow(d))
	after_DCCAS_evals = which(d$obsdate_2m > dccas_end_date)
	nothreat_m2[-after_DCCAS_evals] = bizdays(d$obsdate_2m[-after_DCCAS_evals],dccas_end_date, cal)
	nothreat_m2_all = bizdays(d$obsdate_2m,rep(eval_end_date,nrow(d)),cal)
	nothreat_m2[after_DCCAS_evals] = 0

	# Total nothreat_me
	nothreat_total_m = nothreat_m1+nothreat_m2

	max_ntm2_tot = rep(bizdays(cutoff_m1, eval_end_date, cal), nrow(d))

	# nothreat ME1 that occurred in window 1 -> only until the start of P2 window
	ntm1_cycle1 = rep(NA, nrow(d))
	ntm1_cycle1[-late_tests] = bizdays(d$obsdate_1m[-late_tests], rep(cutoff_p1,nrow(d))[-late_tests], cal)
	ntm1_cycle1[which(ntm1_cycle1 < 0)] = 0

	# nothreat ME2 that occurred in window 4 -> only until the end of p2 window
	ntm2_cycle4 = rep(NA, nrow(d))
	ntm2_cycle3 = rep(NA, nrow(d))
	ntm2_cycle3[-after_DCCAS_evals] = bizdays(d$obsdate_2m[-after_DCCAS_evals], rep(cutoff_p2,nrow(d))[-after_DCCAS_evals], cal)
	ntm2_cycle3[which(ntm2_cycle3 < 0)] = 0
	ntm2_cycle4 = nothreat_m2 - ntm2_cycle3
	
	#################################
	#
	# Admin NO THREAT - 3 CYCLES
	#
	#################################

	# Admin cycle 1
	late_tests = which(d$obsdate_1p > cutoff_p1)
	nothreat_p1 = rep(NA, nrow(d))
	nothreat_p1[-late_tests] = bizdays(d$obsdate_1p[-late_tests], cutoff_p1, cal)

	max_ntp1 = 	rep(bizdays(eval_start, cutoff_p1, cal), nrow(d))

	# Admin cycle 2
	late_tests = which(d$obsdate_2p > cutoff_p2)
	nothreat_p2 = rep(NA, nrow(d))
	nothreat_p2[-late_tests] = bizdays(d$obsdate_2p[-late_tests], cutoff_p2, cal)

	max_ntp2 = 	rep(bizdays(cutoff_p1, cutoff_p2, cal), nrow(d))

	# specify pre and postbreak p2 that occurs within cycle2 only
	max_ntp2c2_postbreak = max_ntm1_postbreak #this is max p2 days postbreak pre cycle3
	max_ntp2c3 = bizdays(cutoff_m1, cutoff_p2, cal) # max nothreat_p2 days in cycle3
	ntp2c2_prebreak = nothreat_p2 - max_ntp2c2_postbreak - max_ntp2c3
	ntp2c2_prebreak[which(ntp2c2_prebreak < 0)] = 0
	ntp2c2_postbreak = nothreat_p2 - max_ntp2c3
	ntp2c2_postbreak[which(ntp2c2_postbreak < 0)] = 0

	# specify pre and post p2 regardless of cycle2
	max_ntp2_postbreak = bizdays(end_winter_break, cutoff_p2, cal) # max nothreat_p2
	nothreat_p2_postbreak = nothreat_p2
	nothreat_p2_postbreak[which(nothreat_p2_postbreak > max_ntp2_postbreak)] = max_ntp2_postbreak
	
	# Admin cycle 3
	nothreat_p3 = rep(NA, nrow(d))
	nothreat_p3_all = rep(NA, nrow(d))

	after_DCCAS_evals = which(d$obsdate_3p > dccas_end_date)
	after_end_of_school = which(d$obsdate_3p > eval_end_date)
	
	if(length(after_end_of_school)>0){
		nothreat_p3_all[-after_end_of_school] = bizdays(d$obsdate_3p[-after_end_of_school], eval_end_date, cal)
		nothreat_p3_all[after_end_of_school] = 0
	} else {
		nothreat_p3_all = bizdays(d$obsdate_3p, rep(eval_end_date,nrow(d)), cal)
	}
	nothreat_p3[-after_DCCAS_evals] = bizdays(d$obsdate_3p[-after_DCCAS_evals],dccas_end_date, cal)
	nothreat_p3[after_DCCAS_evals] = 0


	max_ntp3_tot = rep(bizdays(cutoff_p2, eval_end_date, cal), nrow(d))
	max_ntp3_pretest = rep(bizdays(cutoff_p2, dccas_end_date, cal), nrow(d))

	# Total nothreat_p
	nothreat_total_p = nothreat_p1+nothreat_p2+nothreat_p3


	#################################
	#
	# TOTAL NO THRAT - 4 CYCLES
	#
	#################################
	
	# CYCLE 1
	last_test_date_cycle1 = apply(d[,c("obsdate_1m","obsdate_1p")],1,max)
	nothreat_cycle1 = rep(NA, nrow(d))
	
	positive_no_threat = which(last_test_date_cycle1 <= cutoff_p1)
	not_late_m_eval = which(d$obsdate_1m <= cutoff_m1) # this will also not include those missing me eval
	not_late_p_eval = which(d$obsdate_1p <= cutoff_p1)
	zero_no_threat = intersect(which(last_test_date_cycle1 > cutoff_p1), not_late_m_eval)
	
	nothreat_cycle1[positive_no_threat] = bizdays(last_test_date_cycle1[positive_no_threat],cutoff_p1,cal)
	nothreat_cycle1[zero_no_threat] = 0
	if(this_year == 2010) {
	  nothreat_cycle1[positive_no_threat] = 0
	}


	# CYCLE 2
	if (this_year == 2010){
	  last_test_date_cycle2 = d[,"obsdate_2p"]
	} else {
	  last_test_date_cycle2 = apply(d[,c("obsdate_1m","obsdate_2p")],1,max)
	}
	nothreat_cycle2 = rep(NA, nrow(d))
	
	ntc2_prebreak = rep(NA, nrow(d))
	ntc2_postbreak = rep(NA, nrow(d))

	max_postbreak_days = bizdays(end_winter_break,cutoff_m1,cal)

	positive_no_threat = which(last_test_date_cycle2 <= cutoff_m1)
	not_late_m_eval = which(d$obsdate_1m <= cutoff_m1)
	not_late_p_eval = which(d$obsdate_2p <= cutoff_p2)
	zero_no_threat = intersect(intersect(which(last_test_date_cycle2 > cutoff_m1), not_late_m_eval), not_late_p_eval)
	
	nothreat_cycle2[positive_no_threat] = bizdays(last_test_date_cycle2[positive_no_threat],cutoff_m1,cal)
	nothreat_cycle2[zero_no_threat] = 0

	ntc2_prebreak[zero_no_threat] = 0
	ntc2_postbreak[zero_no_threat] = 0

	ntc2_prebreak[positive_no_threat] = nothreat_cycle2[positive_no_threat] - max_postbreak_days
	ntc2_prebreak[which(ntc2_prebreak < 0)] = 0
	ntc2_postbreak[positive_no_threat] = nothreat_cycle2[positive_no_threat] - ntc2_prebreak[positive_no_threat]
	
	max_ntc2 = rep(bizdays(cutoff_p1, cutoff_m1, cal), nrow(d))

	
	
	# CYCLE 3
	last_test_date_cycle3 = apply(d[,c("obsdate_2m","obsdate_2p")],1,max)
	nothreat_cycle3 = rep(NA, nrow(d))
	
	positive_no_threat = which(last_test_date_cycle3 <= cutoff_p2)
	not_late_m_eval = which(d$obsdate_2m <= eval_end_date)
	not_late_p_eval = which(d$obsdate_2p <= cutoff_p2)
	zero_no_threat = intersect(intersect(which(last_test_date_cycle3 > cutoff_p2), not_late_m_eval),not_late_p_eval)
	
	nothreat_cycle3[positive_no_threat] = bizdays(last_test_date_cycle3[positive_no_threat],cutoff_p2,cal)
	nothreat_cycle3[zero_no_threat] = 0

	max_ntc3 = rep(bizdays(cutoff_m1, cutoff_p2, cal), nrow(d))
	
	
	
	# CYCLE 4
	last_test_date_cycle4 = apply(d[,c("obsdate_2m","obsdate_3p")],1,max)
	after_DCCAS_evals = which(last_test_date_cycle4 > dccas_end_date)
	not_late_final_eval = which(last_test_date_cycle4 <= eval_end_date)
	after_dccas_but_not_late = intersect(not_late_final_eval, after_DCCAS_evals)
	
	nothreat_cycle4 = rep(NA, nrow(d))
	nothreat_cycle4[-after_DCCAS_evals] = bizdays(last_test_date_cycle4[-after_DCCAS_evals],dccas_end_date)
	nothreat_cycle4[after_dccas_but_not_late] = rep(0,length(after_dccas_but_not_late))
	
	# No-Threat Cycle 4 ignore test
	not_late_final_eval = which(last_test_date_cycle4 <= eval_end_date)
	
	posttest_ntc4 = rep(NA, nrow(d))
	posttest_ntc4[not_late_final_eval] = bizdays(last_test_date_cycle4[not_late_final_eval],eval_end_date)
	posttest_ntc4[-not_late_final_eval] = 0
	
	
	#### LATE VARIABLES ###
	
	late_p1 = d[,"obsdate_1p"] > cutoff_p1
	late_p2 = d[,"obsdate_2p"] > cutoff_p2
	late_p3 = d[,"obsdate_3p"] > eval_end_date
	
	late_m1 = d[,"obsdate_1m"] > cutoff_m1
	late_m2 = d[,"obsdate_2m"] > eval_end_date
	
	

	
	
	
	# TOTAL NO THREAT
	nothreat_total = nothreat_cycle1 + nothreat_cycle2 + nothreat_cycle3 + nothreat_cycle4

	#################################
	#
	# post-conference time till test
	#
	#################################
	postconf = as.data.frame(matrix(NA, ncol=0,nrow=nrow(d)))

	postconf$postconf_m1 = bizdays(d[,"confdate_1m"], rep(dccas_start_date, nrow(d)), cal)
	postconf$postconf_m2 = bizdays(d[,"confdate_2m"], rep(dccas_start_date, nrow(d)), cal)
	postconf$postconf_p1 = bizdays(d[,"confdate_1p"], rep(dccas_start_date, nrow(d)), cal)
	postconf$postconf_p2 = bizdays(d[,"confdate_2p"], rep(dccas_start_date, nrow(d)), cal)
	postconf$postconf_p3 = bizdays(d[,"confdate_3p"], rep(dccas_start_date, nrow(d)), cal)

	# replace postconf days that are negative with 0
	postconf$postconf_m2[which(postconf$postconf_m2 < 0)] = 0
	postconf$postconf_p3[which(postconf$postconf_p3 < 0)] = 0

	# Break up m1 and p2 into pre- and post-break post-conference time - is there a reason to do this?
#	max_m1_days_postbreak = bizdays(end_winter_break,dccas_start_date, cal)
#	postconf$postconf_m1_prebreak = postconf$postconf_m1 - max_m1_days_postbreak
#	postconf$postconf_m1_prebreak[which(postconf$postconf_m1_prebreak < 0)] = 0
#	postconf$postconf_m1_postbreak = postconf$postconf_m1 - postconf$postconf_m1_prebreak

#	max_p2_days_postbreak = bizdays(end_winter_break,dccas_start_date,cal)
#	postconf$postconf_p2_prebreak = postconf$postconf_p2 - max_p2_days_postbreak
#	postconf$postconf_p2_prebreak[which(postconf$postconf_p2_prebreak < 0)] = 0
#	postconf$postconf_p2_postbreak = postconf$postconf_p2 - postconf$postconf_p2_prebreak

	#################################
	#
	# Evaluation probability data
	#
	#################################

	# For each school, we need the daily probability of ME and Principal evaluation
	# For each observation, we need the probability of evaluation on that day
	# For each teacher, we need the average professional accountability for each potential observation
	
	teacher_school = rep(NA, nrow(d))
	for(t in 1:length(teacher_school))
	{
		s = as.character(d2$sch1[which(d2$id == d$id[t])])
		if(length(s)>0)
		{
			teacher_school[t] = s
		}
	}

	total.evals.school = table(teacher_school)
	all_schools = names(total.evals.school)
	days.m1 = bizseq(academic_calendar$important[which(academic_calendar$name=="Evaluation Start Date")],
				academic_calendar$important[which(academic_calendar$name=="Master Educator Cutoff 1")], cal)

	days.m2 = bizseq(academic_calendar$important[which(academic_calendar$name=="Master Educator Cutoff 1")],
				academic_calendar$important[which(academic_calendar$name=="Evaluation End Date")], cal)

	days.p1 = bizseq(academic_calendar$important[which(academic_calendar$name=="Evaluation Start Date")],
				academic_calendar$important[which(academic_calendar$name=="Administrative Cutoff 1")], cal)

	days.p2 = bizseq(academic_calendar$important[which(academic_calendar$name=="Administrative Cutoff 1")],
				academic_calendar$important[which(academic_calendar$name=="Administrative Cutoff 2")], cal)

	days.p3 = bizseq(academic_calendar$important[which(academic_calendar$name=="Administrative Cutoff 2")],
				academic_calendar$important[which(academic_calendar$name=="Evaluation End Date")], cal)

	eval.counts.m1 = matrix(NA, nrow=length(total.evals.school), ncol=length(days.m1), dimnames=list(all_schools))
	eval.counts.m2 = matrix(NA, nrow=length(total.evals.school), ncol=length(days.m2), dimnames=list(all_schools))
	eval.counts.p1 = matrix(NA, nrow=length(total.evals.school), ncol=length(days.p1), dimnames=list(all_schools))
	eval.counts.p2 = matrix(NA, nrow=length(total.evals.school), ncol=length(days.p2), dimnames=list(all_schools))
	eval.counts.p3 = matrix(NA, nrow=length(total.evals.school), ncol=length(days.p3), dimnames=list(all_schools))

	all.evalcounts = list(eval.counts.m1, eval.counts.m2, eval.counts.p1, eval.counts.p2, eval.counts.p3)
	all.evalremaining = list(eval.counts.m1, eval.counts.m2, eval.counts.p1, eval.counts.p2, eval.counts.p3)
	all.cumevals = list(eval.counts.m1, eval.counts.m2, eval.counts.p1, eval.counts.p2, eval.counts.p3)
	all.days = list(days.m1, days.m2, days.p1, days.p2, days.p3)
	columns = c("obsdate_1m","obsdate_2m","obsdate_1p","obsdate_2p","obsdate_3p")

	for(s in all_schools)
	{	
		teacher.ids = d2$id[which(d2$sch1 == s)]
		this.school = d[which(d$id %in% teacher.ids),]
		for(m in 1:length(all.evalcounts))
		{
			this.days = all.days[[m]]
			obs.col = columns[m]
			for(t in 1:length(this.days))
			{
				this.date = this.days[t]
				evals.total = nrow(this.school)
				evals.complete = sum(this.school[,obs.col]<this.date, na.rm=TRUE)
				evals.remaining = sum(this.school[,obs.col]>=this.date, na.rm=TRUE)
				evals.today = sum(this.school[,obs.col]==this.date, na.rm=TRUE)

				all.evalcounts[[m]][s,t] = evals.today
				all.evalremaining[[m]][s,t] = evals.remaining
				all.cumevals[[m]][s,t] = evals.complete
			}
		}
	}

	## Calculate all cumulative hazard values
	hazard.list = list()
	hazard.uniform = list()
	prob.list = list()
	for(m in 1:5)
	{
		inner_hazard = matrix(NA, ncol=ncol(all.evalcounts[[m]]), nrow=nrow(all.evalcounts[[m]]))
		rownames(inner_hazard) = all_schools

		inner_hazard.uniform = inner_hazard
		inner_prob = inner_hazard

		for(s in all_schools)
		{
			results = calculate.hazard(m=m, s=s, all.evalcounts=all.evalcounts, all.evalremaining=all.evalremaining,
				all.cumevals=all.cumevals, columns=columns, type="kernel", plot=FALSE)
			results.uniform = calculate.hazard(m=m, s=s, all.evalcounts=all.evalcounts, all.evalremaining=all.evalremaining,
				all.cumevals=all.cumevals, columns=columns, type="uniform", plot=FALSE)

			if(!is.null(results))
			{
				# take the cumulative probabilities up to each day (CDF)
				# then when applied to teacher, their total prob is the CDF of their obs date
				inner_hazard[s,] = results$sum.p
				inner_hazard.uniform[s,] = results.uniform$sum.p
				inner_prob[s,] = results$daily.p
			}
		}
	
		hazard.list[[length(hazard.list)+1]] = inner_hazard
		hazard.uniform[[length(hazard.uniform)+1]] = inner_hazard.uniform
		prob.list[[length(prob.list)+1]] = inner_prob

		names(hazard.list)[length(hazard.list)]=columns[m]
		names(hazard.uniform)[length(hazard.uniform)]=columns[m]
		names(prob.list)[length(prob.list)]=columns[m]
	}

	## assign cumulative hazard values to teachers
	## for each teacher record, there is a hazard value for each observation (5 in all)
	## so need txo identify which day of period evaluation occurred because this will 
	## identify the column in the hazard matrix.

	all.teacher.haz = list(haz_m1 = rep(NA, nrow(d)), haz_m2 = rep(NA, nrow(d)), 
		haz_p1 = rep(NA, nrow(d)), haz_p2 = rep(NA, nrow(d)), haz_p3 = rep(NA, nrow(d)))
	all.teacher.haz.uniform = list(haz_unif_m1 = rep(NA, nrow(d)), haz_unif_m2 = rep(NA, nrow(d)), 
		haz_unif_p1 = rep(NA, nrow(d)), haz_unif_p2 = rep(NA, nrow(d)), haz_unif_p3 = rep(NA, nrow(d)))

	all.teacher.haz.beyondtest = list(haz_m2_beyond = rep(NA, nrow(d)), haz_p3_beyond = rep(NA, nrow(d)))
	all.teacher.thaz.p = list(haz_m1_totp = rep(NA, nrow(d)), haz_m2_totp = rep(NA, nrow(d)), 
		haz_p1_totp = rep(NA, nrow(d)), haz_p2_totp = rep(NA, nrow(d)), haz_p3_totp = rep(NA, nrow(d)))
	all.teacher.thaz.sum = list(haz_m1_tot_sum = rep(NA, nrow(d)), haz_m2_tot_sum = rep(NA, nrow(d)), 
		haz_p1_tot_sum = rep(NA, nrow(d)), haz_p2_tot_sum = rep(NA, nrow(d)), haz_p3_tot_sum = rep(NA, nrow(d)))


	all.weighted.nothreat = list(ntc1_weighted = rep(NA,nrow(d)), ntc2_weighted = rep(NA,nrow(d)),
		ntc3_weighted = rep(NA,nrow(d)), ntc4_weighted = rep(NA,nrow(d)))

	test = 0
	for(i in 1:nrow(d))
	{
		s = teacher_school[i]
		for(m in 1:5)
		{
			this.days = all.days[[m]]
			obs.col = columns[m]
			day.numeric = which(this.days == d[i,obs.col])
			if(length(day.numeric)>0)
			{
				if(d[i,obs.col]>dccas_end_date){
					day.numeric = which(this.days == dccas_end_date)
				}
				all.teacher.haz[[m]][i] = hazard.list[[m]][s,day.numeric]
				all.teacher.haz.uniform[[m]][i] = hazard.uniform[[m]][s,day.numeric]	
			}
		}

		# create hazard that goes beyond test date
		bridge=list("2"=1,"5"=2)
		for(m in c(2,5))
		{
			index = bridge[[as.character(m)]]
			this.days = all.days[[m]]
			obs.col = columns[m]
			day.numeric = which(this.days == d[i,obs.col])
			if(length(day.numeric)>0)
			{
				max.days = length(hazard.list[[m]][s,])
				if(day.numeric <= max.days){
					all.teacher.haz.beyondtest[[index]][i] = hazard.list[[m]][s,day.numeric]
				} else{
					print(d[i,obs.col])
					test = test + 1
				}
			}
		}

		# create "total hazard" on day of evaluation
		
		## THAZ FUNCTION ##
		get.thaz = function(m1, m2, m3, prob.data){
			# return NA if eval was outside range
			c1 = columns[m1]
			if(!(d[i,c1] %in% all.days[[m1]]) | is.na(d[i,c1])){return(NA)}

			start.dates = list("1"=eval_start,"2"=cutoff_m1,"3"=eval_start,"4"=cutoff_p1,"5"=cutoff_p2)
			end.dates = list("1"=cutoff_m1,"2"=eval_end_date,"3"=cutoff_p1,"4"=cutoff_p2,"5"=eval_end_date)
			
			c2 = columns[m2]
			m2.start = start.dates[[as.character(m2)]]
			miss2 = is.na(d[i,c2])

			if(!is.null(m3)){
				c3 = columns[m3]
				m3.start = start.dates[[as.character(m3)]]
				miss3 = is.na(d[i,c3])
			}
			
			h1 = prob.data[[m1]][s,which(all.days[[m1]] == d[i,c1])]
			h2 = 0
	
			if(!is.null(m3)){
				if(!miss2 && d[i,c1] < d[i,c2] && d[i,c1] > m2.start){
					if(d[i,c1] %in% all.days[[m2]]){
						h2 = prob.data[[m2]][s,which(all.days[[m2]] == d[i,c1])]
					}
				} else if(!miss3 && d[i,c1] < d[i,c3] && d[i,c1] > m3.start){
					if(d[i,c1] %in% all.days[[m3]]){
						h2 = prob.data[[m3]][s,which(all.days[[m3]] == d[i,c1])]
					}
				}
			} else {
				if(!miss2 && d[i,c1] < d[i,c2] && d[i,c1] > m2.start){
					if(d[i,c1] %in% all.days[[m2]]){
						h2 = prob.data[[m2]][s,which(all.days[[m2]] == d[i,c1])]
					}
				}
			}
			
			return(h1+h2-h1*h2)
		}

		# ME1 (m = 1) p1 is possible (m=3) and p2 is possible (m=4)
		all.teacher.thaz.p[[1]][i] = get.thaz(m1=1, m2=3, m3=4, prob.list)
		all.teacher.thaz.sum[[1]][i] = get.thaz(m1=1, m2=3, m3=4, hazard.list)

		# ME2 (m = 2) p2 is possible (m=4) and p3 is possible (m=5)
		all.teacher.thaz.p[[2]][i] = get.thaz(m1=2, m2=4, m3=5, prob.list)
		all.teacher.thaz.sum[[2]][i] = get.thaz(m1=2, m2=4, m3=5, hazard.list)

		# P1 (m = 3) m1 is possible (m=1)
		all.teacher.thaz.p[[3]][i] = get.thaz(m1=3, m2=1, m3=NULL, prob.list)
		all.teacher.thaz.sum[[3]][i] = get.thaz(m1=3, m2=1, m3=NULL, hazard.list)

		# P2 (m = 4) m1 is possible (m=1) and m2 is possible (m=2)
		all.teacher.thaz.p[[4]][i] = get.thaz(m1=4, m2=1, m3=2, prob.list)
		all.teacher.thaz.sum[[4]][i] = get.thaz(m1=4, m2=1, m3=2, hazard.list)

		# P3 (m = 5) m2 is possible (m=2)
		all.teacher.thaz.p[[5]][i] = get.thaz(m1=5, m2=2, m3=NULL, prob.list)
		all.teacher.thaz.sum[[5]][i] = get.thaz(m1=5, m2=2, m3=NULL, hazard.list)

		
	}



		# function to calculate weighted nothreat
############ COMMENTED OUT  ##################
		if(2>1){
		get_weighted_nothreat = function(limiting.column,other.column,end_date,person,nothreat_data,s,columns)
		{
			
			if(!is.na(nothreat_data))
			{
				if(nothreat_data > 0)
				{
					start.no.threat.limiting = which(all.days[[limiting.column]] == 
						max(d[person,columns[limiting.column]], d[person,columns[other.column]]))
					start.no.threat.other = which(all.days[[other.column]] == 
						max(d[person,columns[limiting.column]], d[person,columns[other.column]]))
					end.no.threat.other = which(all.days[[other.column]] == end_date)
	
					start.hazard.limiting = hazard.list[[limiting.column]][s,start.no.threat.limiting]
					end.hazard.limiting = max(hazard.list[[limiting.column]],na.rm=TRUE)
					
					start.hazard.other = hazard.list[[other.column]][s,start.no.threat.other]
					end.hazard.other = hazard.list[[other.column]][s,end.no.threat.other]
	
					return_value = end.hazard.limiting - start.hazard.limiting + 
						end.hazard.other - start.hazard.other
				} else {return_value = 0}
			} else {return_value = NA}
			return(return_value)
		}
		}

		# p1 limits (3), m1 is other (1)
		#all.weighted.nothreat$ntc1_weighted[i] = get_weighted_nothreat(3,1,cutoff_p1,i,nothreat_cycle1[i],s,columns)

		# m1 limits (1), p2 is other (4)
		#all.weighted.nothreat$ntc2_weighted[i] = get_weighted_nothreat(1,4,i,cutoff_m1,nothreat_cycle2[i],s,columns)

		# p2 limits (4), m2 is other (2)
		#all.weighted.nothreat$ntc3_weighted[i] = get_weighted_nothreat(4,2,i,cutoff_p2,nothreat_cycle3[i],s,columns)

		# m2 limits (2), p3 is other (5)
		#all.weighted.nothreat$ntc4_weighted[i] = get_weighted_nothreat(2,5,i,dccas_enddate,nothreat_cycle4[i],s,columns)
######################################

	


	all.teacher.haz = as.data.frame(all.teacher.haz)
	all.teacher.haz.uniform = as.data.frame(all.teacher.haz.uniform)


	#par(mfrow=c(5,1),mar=c(2,5,2,2))
	#for(m in 1:5)
	#{
		#plot(density(na.omit(all.teacher.haz[,m])))
	#}

	#dev.new()
	#par(mfrow=c(5,1),mar=c(2,5,2,2))
	#for(m in 1:5)
	#{
		#plot(density(na.omit(all.teacher.haz[,m])))
	#}



	
	#################################
	#
	# Save Data
	#
	#################################
	
	threat.data = cbind(as.numeric(d$id), as.numeric(d$year), as.data.frame(nothreat_total_m), 
		as.data.frame(nothreat_m1), as.data.frame(ntm1_prebreak), as.data.frame(ntm1_postbreak), 
		as.data.frame(ntm1_cycle1), as.data.frame(nothreat_m2), as.data.frame(nothreat_m2_all), as.data.frame(nothreat_p3_all),
		as.data.frame(nothreat_total_p), as.data.frame(nothreat_p1), as.data.frame(nothreat_p2), as.data.frame(nothreat_p2_postbreak), as.data.frame(nothreat_p3), as.data.frame(max_ntp3_pretest), 
		as.data.frame(nothreat_cycle1),as.data.frame(nothreat_cycle2), as.data.frame(max_ntc2), as.data.frame(max_ntc3),
		as.data.frame(ntc2_prebreak),as.data.frame(ntc2_postbreak),	as.data.frame(nothreat_cycle3),
		as.data.frame(ntp2c2_prebreak),as.data.frame(ntp2c2_postbreak), as.data.frame(ntm2_cycle3), as.data.frame(ntm2_cycle4),
		as.data.frame(max_ntp1),as.data.frame(max_ntp2),as.data.frame(max_ntp3_tot),as.data.frame(max_ntm1),as.data.frame(max_ntm2_tot),
		as.data.frame(nothreat_cycle4),as.data.frame(nothreat_total),as.data.frame(posttest_ntc4) ,
		as.data.frame(late_p1),as.data.frame(late_p2),as.data.frame(late_p3),as.data.frame(late_m1),as.data.frame(late_m2),
		all.teacher.haz, all.teacher.haz.uniform, postconf, all.teacher.haz.beyondtest,
		  all.teacher.thaz.p, all.teacher.thaz.sum)
	names(threat.data)[c(1,2)] = c("id","year")
	write.dta(threat.data, paste(local_data_directory,"\\build\\temp\\threat_", this_year, ".dta",sep=''))

	#export.data = convert_string_to_factor(cbind(d,as.data.frame(threat.data)))
	#write.dta(export.data, paste(local_data_directory,"build\\threat_", this_year, ".dta",sep=''))
	#write_dta(export.data, paste(local_data_directory,"build\\temp\\threat_test", this_year, ".dta",sep=''))
}


